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INTRODUCTION 


A  common  occurrence  in  battlefield  situations  is  the  loss  of  meteorologi¬ 
cal  (met)  data  due  to  malfunctions  or  loss  of  weather  balloons.  Recourse  must 
be  made  to  other  existing  balloon  data,  both  contemporary  and  dated.  For 
artillery  purposes  a  complete  met  message  is  required.  In  this  report  we 
describe  a  least  squares  regressional  analysis  approach  to  supplying  the 
missing  meteorological  data  necessary  to  complete  the  met  message. 

Using  data  from  the  PASS  (Prototype  Artillery  Subsystem)  field  experiment 
(1)  conducted  at  White  Sands  Missile  Range,  New  Mexico,  during  October  - 
December,  1974,  we  have  formulated  several  least  squares  routines  to  predict 
four  variables  from  available  data  -  wind  velocity  and  direction,  temperature, 
and  pressure.  An  artillery  shot  was  simulated  on  a  computer  by  calculating  an 
appropriate  ballistic  trajectory  from  the  actual  meteorological  data  for 
layers  0-4  from  one  balloon  combined  with  predicted  values  of  missing  data  'or 
layers  5-9  supplied  by  each  test  routine.  These  results  were  compared  to  the 
results  obtained  from  using  the  entire  set  of  actual  meteorological  data  as 
recorded  by  the  balloon.  In  turn,  each  of  the  methods  were  intercompared  to 
identify  the  best  candidates  for  further  comparisons  with  methods  currently  in 
use. 


THE  LEAST  SQUARES  METHOD 

Each  simulated  shot  was  aimed  at  a  target  9500  meters  north  of  the  bat¬ 
tery,  with  a  trajectory  apex  of  4000  meters  (layer  9).  A  155  mm  howitzer  with 
a  7W  charge  was  used  in  each  case,  located  at  the  met  site.  Using  met  data 
from  surrounding  sites  and  from  layers  0-4  at  the  local  site,  a  prediction  was 
made  for  the  meteorological  parameters  at  layers  5-9.  Following  procedures 
from  the  firing  tables  (2)  using  artillery  met  messages  (3),  the  artillery 
trajectory  was  calculated  and  compared  to  one  made  with  the  actual  data  from 
layers  0-9  as  measured  by  the  local  balloon. 

Four  met  sites  from  the  PASS  experiment  were  used.  Their  names,  coordi¬ 
nates,  and  elevations  are  listed  in  Table  1.  The  artillery  firings  were  made 
only  from  SMR  and  MCG.  In  order  to  sample  conditions  throughout  the  day, 
firings  were  made  every  hour  (±  15  mins)  beginning  at  0600  (local  time)  and 
ending  at  1600.  Table  2  shows  the  dates  and  times  of  the  firings.  For  each 
simulated  firing  it  was  necessary  to  have  complete  data  from  each  of  the  four 


met  stations  for  both  the  time  of  the  firing  and  two  hours  previous.  This 
severely  restricted  the  data  available  for  analysis,  which  resulted  in  having 
to  select  each  two  hour  block  of  data  from  a  different  day.  For  this  reason  a 
detailed  direct  comparison  of  results  as  a  function  of  time  of  day  is  not 
possible,  but  because  the  weather  conditions  were  similar  and  stable  through¬ 
out  the  period,  we  do  make  some  inferences  with  regard  to  time  of  day. 


Table  1.  White  Sands  transverse  mercator  system  coordinates 
of  met  stations  (converted  to  meters) 

X  (increases  Y  (increases  Z  (above  mean 
Station  Name _ to  the  East)  to  the  North) _ sea  level) 


Small  Missile  Range 

(SMR) 

144040 

65614 

1219 

McGregor  (MCG) 

165731 

42786 

1249 

Orogrande  (0R0) 

169896 

57769 

1280 

LC-36  (LSX) 

153348 

57831 

1229 

Table  2. 

Date  and  Time 

of  Selected  Data 

Date  1974 

Local  Stan 

11/14 

0545 

11/11 

0645 

11/12 

0815 

11/11 

0845 

11/15 

1000 

12/02 

1115 

11/14 

1145 

11/27 

1315 

11/20 

1345 

11/08 

1445 

11/20 

1545 

For  each  simulated  trajectory  computed  with  actual/predicted  data  we  con¬ 
structed  a  measure  of  error  attributable  to  the  cross  wind  component  (V^), 
head  wind  component  (V  ) ,  temperature  (t),  and  density  (D).  (The  balloon- 
measured  pressure  in  millibars  and  temperature  in  degrees  Kelvin  are  used  to 
find  the  density:  D  =  348.4  P/t).  This  measure  (A)  is  the  difference  in 
corrections  determined  from  the  firing  tables  converted  to  meters  between  the 
method  using  actual/predicted  data  and  the  control  method  using  actual  data 
for  all  ten  levels.  Combining  these  errors  we  then  express  a  bias  (0)  and 
variance  (a2)  for  each  method.  The  total  miss  distance  or  bias  is  the  vector 
sum  of  the  cross  and  range  miss  distances: 


6 


p  =  [AV2  +  (AV  +  At  +  AD)2]5*, 
x  y 

The  total  error  squared  or  variance  is  the  sum  of  the  square  of  the  components 

a2  =  AV2  +  AV2  +  At2  +  AD2, 
x  y 

Eleven  different  methods  were  tested.  In  each  method  four  separate  least 

squares  analyses  of  the  data  were  made,  one  for  each  variable  V  ,  V  ,  t,  and 

x  y 

In  P,  the  latter  two  variables  then  being  combined  to  form  the  density. 
Differences  (A)  between  computed  and  actual  values  were  combined  to  form  p  and 
or2 . 

In  general  we  can  express  the  fitting  equation  for  a  variable  (say  tem¬ 
perature)  as  a  function  of  position  (X  and  Y) ,  altitude  (Z),  and  time  (T): 


NS  NZ  NT 

(1)  t  =  a  +  .2  (b  X1  +  c.Y1)  +  ,I.d.ZJ  +  JLf.T*. 

i=l  i  l  j=l  j  k-1  k 

If  a  limit  is  zero  (i.e.,  NS,  NT,  or  NZ  is  zero),  then  that  variable  is 
not  included  in  the  analysis.  Vertical  fits  emphasize  higher  order  terms  in  Z 
(NZ>1;  NZ>NS  and  NT)  and  horizontal  fits  allow  higher  orders  in  X  and  Y  (NS>1; 
NZ=0).  Higher  order  terms  in  time  were  not  allowed  because  only  two  times 
were  used  in  the  observational  equations  of  condition,  Tq  and  Tq  -  2  (hrs). 
Since  only  four  stations  were  used  at  most,  NS  was  never  more  than  3,  and  the 
maximum  order  of  Z  was  never  greater  than  5  (NZ<5).  Although  there  were  often 
enough  equations  of  condition  to  accomodate  higher  terms  in  Z  for  vertical 
fits,  preliminary  runs  indicated  that  little  was  gained  in  expanding  the  order 
of  Z  beyond  5.  In  addition,  the  vertical  fit  for  In  P  was  made  only  over  Z, 
and  not  over  X,  Y,  T,  or  higher  powers  of  Z.  Preliminary  attempts  to  include 
them  in  the  fitting  of  In  P  (or  P)  did  not  indicate  an  improvement  over  ignor¬ 
ing  them.  Consequently,  for  vertical  fits  of  In  P,  NS  =  NT  =  0,  and  NZ  =  1, 
in  all  cases. 

For  vertical  fits  to  a  parameter,  one  equation  with  properly  determined 
coefficients  (a-f^)  suffices  to  yield  the  five  missing  values  of  the  variable. 
For  instance,  once  the  least  squares  fitting  equation  (1)  is  solved  for  the 
unknown  coefficients,  the  temperature  can  be  predicted  for  layers  5-9.  For 
horizontal  fits,  however,  each  fitting  equation  must  be  solved  at  each  layer. 
Thus  it  takes  five  such  equations  to  determine  five  missing  temperatures. 
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Many  more  equations  of  conditions  are  available  for  a  vertical  fit  (IS  for  two 
stations)  than  for  a  horizontal  fit  (3  equations  of  condition  for  each  layer). 

Thus  a  two  station  vertical  fit  can  accomodate  15  coefficients,  while  a  hori¬ 
zontal  fit  can  only  accomodate  3  per  layer,  and  must  be  determined  separately 
for  each  of  the  layers  where  data  is  missing.  The  coefficients  in  a  vertical 
fit  are  over  determined,  while  in  many  cases  one  must  choose  which  of  the 
parameters  X,  Y,  or  T  must  be  left  out  in  a  horizontal  fit.  Nevertheless, 
since  many  meteorological  conditions  are  partitioned  into  layers  as  manifest 
by  wind  shears  and  inversions,  a  horizontal  fit  is  appropriate  in  some  cases. 

An  important  characteristic  of  our  approach  is  that  we  take  the  position 
of  the  balloon  into  account  when  expressing  values  of  meteorological  param¬ 
eters.  The  value  of  the  variable  In  P,  for  example,  is  recorded  at  the 
balloon's  actual  position,  assuming  a  rise  rate  of  300  meters/min  and  using 
the  observed  winds  to  track  the  balloon  from  one  layer  to  the  next.  The 
missing  data  is  supplied  by  the  fitting  equation  for  the  position  of  the 
balloon,  which  facilitates  comparison  of  the  actual  to  predicted  data  at  the 
same  position  and  altitude. 

The  eleven  methods  tested  involved  one,  two,  or  four  met  stations.  Each 
station  had  to  have  a  complete  set  of  data  for  two  hours  previous  to  and  for 
the  time  of  the  firing.  The  missing  data  was  simulated  by  simply  ignoring  the 
upper  5  layers  of  the  station  containing  the  fictitious  battery.  The  values 
of  NS,  NT,  and  NZ  in  equation  (1)  are  shown  below  for  each  method  a  through  k. 

1  station 

a)  Vertical:  NS  =  0,  NT  =  1,  NZ  =  5.  (NS=NT=0,  NZ=1,  for  In  P) 

b)  Tacfire:  Upper  5  levels  at  Tq  are  assumed  to  be  the  same  as  upper  levels  at 

Tq-2,  but  adjusted  by  the  difference  between  variables  at  level  4. 

c)  Persistence:  Upper  5  levels  at  Tq  are  assumed  to  be  the  same  as  upper  levels 

at  Tq-2,  without  any  adjustments. 

2  stations 

d)  Vertical:  NS  =  1 ,  NT  =  1 ,  NZ  =  5 .  (NS=NT=0,  NZ  =  1,  for  In  P) . 

e)  Horizontal:  NS  =  0,  NT  =  1 ,  NZ  =  0. 

f)  Horizontal:  NS  =  1,  NT  =  0,  NZ  =  0. 

g)  Vertical  and  Horizontal:  Same  as  d  for  temperature  and  In  P; 

same  as  e  for  V  and  V  . 

X  y 

h)  Vertical  and  Horizontal:  Same  as  d  for  temperature  and  In  P; 

same  as  f  for  winds. 
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4  stations 

i)  Vertical:  NS  =  3,  NT  =  1,  NZ  =  5.  (NS=NT=0,  NZ=1,  for  In  P). 

j)  Horizontal:  NS  =  2,  NT  -  1,  NZ  =  0. 

k)  Vertical  and  Horizontal:  Same  as  i  for  temperature  and  In  P; 

same  as  j  for  winds. 

RESULTS 

The  results  of  the  firings  are  shown  in  Figures  1-11,  which  are  graphs  of 
the  bias  (p)  and  standard  error  (o) ,  averaged  between  SMR  and  MCG,  plotted 
against  time  of  day  for  each  method.  Inspection  of  these  figures  reveals  that 
all  methods  experienced  worse  results  during  the  morning  transition  period 
than  later  during  the  day.  The  stability  of  meteorological  conditions  after 
the  transition  period  and  after  the  dissipation  of  surface  based  temperature 
inversions  leads  to  better  predictions  of  missing  meteorological  parameters  by 
all  methods. 

Table  3  is  a  compilation  of  the  average  bias  and  standard  error,  broken 
down  by  method.  Two  numbers  are  shown  for  each  bias  and  standard  error.  The 
smaller  number  is  the  average  recomputed  after  eliminating  values  in  excess  of 
two  standard  deviations  from  the  original  (larger)  averages.  This  is  done  in 
order  to  make  a  fair  comparison  with  methods  f  and  h,  where  attempts  to  use  a 
horizontal  fit  resulted  in  wild  values  for  the  variables  because  only  3  equa¬ 
tions  of  condition  were  used  to  determine  3  unknown  coefficients.  This  will 
be  discussed  further  below. 

Figure  12  depicts  the  two  values  each  of  the  bias  and  standard  error 
before  and  after  elimination  of  those  values  greater  than  two  standard 
deviations.  Because  the  firings  were  made  due  north  with  the  prevailing  winds 
out  of  the  west,  it  was  expected  that  the  results  from  the  SMR  station  should 
have  been  significantly  poorer  than  from  the  MCG  station.  In  fact,  however, 
SMR  outscored  MCG  (lower  bias)  by  1/3,  which  is  interpreted  as  due  to  terrain 
effects.  Therefore,  to  minimize  these  effects  and  place  the  emphasis  on 
techniques  instead,  the  averages  between  MCG  and  SMR  are  ued  to  make  Table  3 
and  Figure  12. 
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Table  3.  Average  bias  and  standard  errors,  (in  meters) 
before  and  after  elimination  of  extreme  values 


Method 

Bias : 

after  (before  elimination) 

Standard  error: 
after  (before) 

a 

37  (40) 

39  (41) 

b 

1  station 

50  (56) 

49  (55) 

c 

31  (33) 

32  (35) 

d 

31  (34) 

31  (34) 

e 

38  (44) 

37  (43) 

f 

2  stations 

44  (128) 

57  (133) 
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35  (38) 

37  (40) 

h 

53  (129) 

55  (131) 

i 

29  (32) 

30  (33) 

j 

4  stations 

25  (28) 

30  (33) 

k 

25  (27) 

30  (31) 

Figures  1-11.  Figures  1  through  11,  corresponding  to  methods  a  through 
k,  are  presented  from  left  to  right,  top  to  bottom  on  the  next  two  pages.  In 
each  figure  the  dotted  line  represents  the  standard  error  (o)  and  the  solid 
line  the  bias  (p)  for  each  method  plotted  against  the  time  of  day.  Note  the 
better  results  after  mid-morning  in  each  case. 


Figure  12.  The  lower  right  figure  on  the  second  page  of  figures  repre¬ 
sents  the  overall  average  standard  error  and  bias  as  a  function  of  method 
before  and  after  elimination  of  the  values  greater  than  two  standard 
deviations.  Here  note  that  methods  f  and  h  occasionally  yield  wild  results. 
This  phenomenon  is  explained  in  the  text. 
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For  a  consideration  of  the  contribution  of  each  of  the  components  V^,  V  , 
temperature,  and  density,  to  the  total  variance  o2,  each  difference  A  was  normal¬ 
ized  by  squaring  it  and  dividing  it  by  the  total  variance.  From  a  total  of  242 
cases  (eleven  methods,  at  eleven  times,  from  two  stations)  the  following  breakdown 
is  presented. 

V  =  49.6%  of  total  variance,  ranging  from  17.6%  on  11/15  at  1000  at  SMR 
y  to  87.3%  on  11/20  at  1345  at  MCG. 

V  =  23.8%  of  total  variance,  ranging  from  4.7%  on  li/20  at  1345  to  MCG 

to  72.2%  on  11/15  at  1000  at  MCG. 

temperature  =  1.6%  of  total  variance,  ranging  from  0.1%  on  many  occasions 
to  7.0%  on  11/27  at  1315  at  SMR. 

Density  =  25.0%  of  total  variance,  ranging  from  2.3%  on  11/11  at  0845  at 
MCG  to  55.3%  on  11/12  at  0815  at  MCG. 


The  contribution  of  each  component  to  the  bias  in  each  of  the  eleven  methods 
can  be  assessed  by  taking  the  average  absolute  value  of  the  difference  A  for  each 
kind  of  horizontal  and  vertical  fit  for  1,  2,  or  4  stations.  Such  a  breakdown 
by  methods  and  type  of  fit  is  presented  in  Table  4.  Parenthetical  values  for 
methods  f  and  h  denote  inclusion  of  wild  values. 


Table  4.  Average  absolute  values  of  A  (meters) 

Method  V  t  D  V 

_ y  _  x 


1  station 


a  (vertical) 

56 

7 

36 

26 

b  (Tacfire) 

79 

5 

37 

37 

c  (persistence) 

50 

4 

22 

27 

stations 
d  (vertical) 

44 

5 

24 

29 

e  (horizontal) 

45 

8 

33 

42 

f  (horizontal) 

100(242) 

12(13) 

49(52) 

38(58) 

g  (horizontal/vertical) 

45 

5 

24 

42 

h  (horizontal/vertical) 

100(242) 

5 

24 

38(58) 

stations 
i  (vertical) 

45 

6 

25 

23 

j  (horizontal) 

45 

5 

30 

27 

k  (horizontal/vertical) 

45 

6 

25 

27 

It  is  from  such  a  breakdown  that  we  selected  a  mixture  of  fits  (horizontal 
for  winds  and  vertical  for  temperature  and  pressure)  to  combine  into  methods 
g,  h,  and  k.  Since  pressure  is  clearly  best  fit  by  the  vertical  methods  we 
decided  to  apply  a  vertical  fit  to  the  temperature  also,  in  order  to  keep  the 
density  in  a  vertical  structure. 

An  analysis  of  the  contribution  of  each  of  the  components  to  the  bias 
indicate  that  except  for  In  P,  each  variable  is  predicted  about  as  well  (or 
better)  with  a  horizontal  fit  as  with  a  vertical  fit,  if  the  number  of  equa¬ 
tions  of  condition  is  greater  (preferably  much  greater)  than  the  number  of 
unknown  coefficients  to  be  found.  Physically,  this  indicates  that  the  winds 
and  temperature  are  stratified  in  the  atmosphere.  Mathematically,  the  poor 
results  for  vertical  fits,  despite  the  greater  number  of  available  equations 
of  condition,  is  a  product  of  the  discontinuities  associated  with  the  inter¬ 
faces  of  the  layers.  In  other  words,  a  least  squares  smooth  fit  breaks  down 
when  a  sometimes  nearly  discontinuous  variable  is  encountered.  The  fact  that 
the  vertical  fits  perform  as  well  as  they  do  can  be  accounted  for  by  the  over 
specification  of  the  unknown  coefficients  in  vertical  fits  alluded  to  earlier. 

DISCUSSION 

Regarding  the  persistence  method  as  a  degenerate  case  of  a  horizontal 
fit,  where  NS  =  NT  =  NZ  =  0,  an  important  finding  of  this  report  is  that 
horizontal  fits  are  to  be  preferred  over  vertical  fits  when  enough  data  exists 
from  enough  stations.  Let  M  be  the  "freedom",  the  difference  between  the 
number  of  equations  of  condition  (N)  and  the  number  of  unknown  coefficients 
(P).  For  each  method,  then,  Table  5  lists  M,  N,  P,  and  Q,  the  number  of  times 
equation  (1)  must  be  solved  to  yield  the  missing  5  values  of  a  given  variable. 

Obviously,  M  is  much  greater  when  employing  a  vertical  fit  than  when 
attempting  a  horizontal  regression,  yet  except  when  M=0,  the  horizontal  is  as 
good  or  better  than  the  vertical  (with  the  exception  of  fitting  In  P  where  the 
vertical  fit  is  clearly  superior). 

Two  conclusions  can  be  male  regarding  the  "best"  method  for  supplying 
missing  data.  First,  the  greater  the  number  of  stations  the  better  are  the 
results.  Provided  that  there  is  enough  freedom  (i.e.,  the  number  of  coeffi¬ 
cients  to  be  found  is  less  than  the  number  of  conditional  equations),  a  greater 
number  of  stations  allows  better  predictions  of  meteorological  parameters. 

The  second  conclusion,  on  the  other  hand,  is  that  the  difference  between  methods 
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Table  6.  "Freedom"  as  a  function  of  method 


Method 

a 

b 

c 

d 

e 

f 

g 

h 

i 

j 

k 


M  =  N 
8  15 

15 
15 

26  35 

1  3 

0  3 

same  as  d 
same  as  d 
62  75 

1  7 

same  as  i 


no  least  squares  fit  is  made 
no  least  squares  fit  is  made 


9  1 

2  5 

3  5 

for  In  P  and  t;  same  as  e  for  winds 

for  In  P  and  t;  same  as  f  for  winds 

13  1 

6  5 

for  In  P  and  t;  same  as  j  for  winds 


and  number  of  stations  is  less  than  might  be  expected.  Of  the  three  single 
station  methods,  persistence  leads  to  the  lowest  average  bias  and  variance. 
In  fact,  it  performs  just  as  well  as  the  best  of  the  two  station  methods,  a 
vertical  fit,  and  only  slightly  poorer  than  the  best  of  the  four  station 
methods . 


SPACE/TIME  VARIABILITY 

Among  the  various  approaches  to  interpolating  or  extrapolating  meteoro¬ 
logical  data  to  predict  missing  data,  a  commonly  considered  technique  involves 
weighting  met  messages  according  to  the  ages  and  distances  of  the  stations. 

Traylor  (la)  examined  four  techniques  for  use  in  the  PASS  experiment,  two  of 
which  (the  "weighted  average"  and  the  "plane  fit")  involved  weighted  messages. 

The  weighted  average  scheme  Was  adopted  and  discussed  further  by  Stenmark  et 
al  (4)  and  Blanco  and  Traylor  (5).  When  weights  are  used,  an  equivalence 
between  space  and  time  variability  must  be  established.  All  of  these  workers 
assumed  that  the  errors  inherent  in  old  or  distant  messages  vary  as  the  square 
root  of  the  time  or  distance  of  the  message,  and  that  the  equivalence  between 
the  two  is  between  12km/hr  and  46km/hr.  Traylor  (la)  and  Blanco  and  Traylor 
(5)  chose  30km/hr  while  Stenmark  et  al  (4)  chose  35km/hr  for  their  scaling  factors 

In  order  to  test  these  assumptions,  with  an  eye  towards  incorporating 
weighting  factors  into  our  least  squares  approach,  we  found  two  times  in  the 
PASS  data  set  when  eight  hours  of  continuous  coverage  was  available  at  each 
station.  For  each  of  eight  stations  (Table  7),  we  used  the  met  message  that 


Table  7.  Separation  between  Coordinates  of  stations 
used  for  time/space  variability  study. 


McGregor 

MCG 

0R0 

Orogrand"' 

15.6 

- 

War  Road 

20.8 

28.7 

LC-36 

19.5 

16.6 

Small  Missile  Range 

31.5 

27.0 

WAR  LSX  SMR  APA  HMS  RAN 

16.5 

21.9  12.1 


Apache 

Holloman  37.8 

Rampart  25.0  39.4 


was  8,  6,  4,  and  2  hours  old  to  make  our  standard  artillery  firing  (i.e., 
target  due  north,  9500m  away,  trajectory  apex  4000m).  We  then  made  a  compari¬ 
son  in  each  case  to  the  results  from  using  the  actual  current  message. 

Figures  13  (November  14,  1200  ±  0015)  and  14  (November  15,  1215  ±  0015) 
are  plots  of  the  average  bias  for  each  of  the  eight  stations  over  time,  with 
the  standard  deviation  from  the  average  illustrated  as  an  error  bar  on  each 
point.  Also  shown  in  each  figure  are  curves  fitted  to  the  first  three  points 
only  (2,  4,  and  6  hours).  Eight  hour  old  data  dropped  toward  lower  biases  com¬ 
pared  to  six  hour  old  data.  Consequently,  eight  hours  were  not  considered  in 
the  fittings.  Table  8  shows  the  numerical  values  of  the  results  of  the  fits. 
The  t2  fit  is  not  shown  in  the  figure. 

Table  8.  Functional  dependence  of  Bias  (B)  over  time 


Function 

Standard  i 

11/14 

18.28  Vt 

1.04 

7.95  t 

13.55 

13.45 

8.46  t 

13.08 

1.60  t2 

25.03 

11/15 

22.84  Vx 

28.53 

8.29  t1-08 

20.62 

10.57  t 

18.81 

2.00  t2 

13.61 
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Figure  13.  -  (Top  of  next  page).  Bias  vs.  time  for  November  14,  1974  at 
1145  for  MCG,  0R0,  WAR,  LSX,  and  SMR,  and  1215  for  APA,  HMS,  and  RAM.  Three 
fits  over  the  first  three  points  are  shown:  a  square  root  fit  is  denoted  by 
the  dot-dash  line,  a  linear  fit  by  the  solid  line,  and  a  power  fit  by  the  dashed 
line.  The  standard  deviation  from  the  mean  is  shown  as  an  error  bar  on  each 
point.  The  total  length  of  the  bar  is  2a.  The  numerical  values  for  the  fits 
are  given  in  Table  7. 

Figure  14.  -  (Bottom  of  next  page).  Same  as  for  Figure  13,  but  for  November 
15,  1974  and  15  minutes  later.  See  Table  8  for  the  numerical  values  for  the  fits. 
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Similarly,  current  met  messages  from  every  other  station  were  used  in 
turn  as  a  substitute  for  the  actual  message  of  each  station.  The  biases  were 
calculated  and  are  plotted  as  a  function  of  station  separation  in  Figures  15 
and  16.  Table  9  is  the  result  of  fitting  all  of  the  distance  points.  The  x2 
fit  again  is  not  shown  in  the  figures  in  order  to  reduce  the  clutter  of  too 
many  lines.  Because  the  HMS,  APA,  and  RAM  sites  released  their  balloons  30 
minutes  after  the  other  five  stations,  they  are  treated  separately.  The  three 
were  not  used  to  supply  a  current  met  message  to  any  of  the  five  stations,  and 
vice-versa . 


Table  i.  Functional  dependence  of  Bias  over  space 


Function 

Standard  i 

11/14 

5.88  -fx 

9.26 

1.09  x  •" 

9.23 

1.13  x 

8.81 

0.037  x2 

13.69 

11/15 

3.64  yfx 

7.56 

1.18  x  ♦80 

8.23 

.70  x 

8.43 

0.023  x2 

11.86 

Examining  Figures  13-16,  it  is  not  easy  to  select  the  best  fit.  The 
residuals  from  the  fits  (expressed  in  Tables  8  and  9  in  the  standard  errors) 
do  not  indicate  a  significant  difference  between  functional  forms  of  variation 
either.  There  is  little  reason  to  choose  a  linear,  square  root,  or  power 
relationship  over  either  distance  or  time.  Only  a  quadratic  fit  can  be  elimi¬ 
nated  as  inferior,  although  it  does  give  the  best  fit  over  time  on  11/15 
according  to  the  standard  errors.  The  assumption  of  a  square  root  dependence 
of  variability  on  time  and  space  does  not  appear  to  be  justified,  but  from  the 
sparse  data  analyzed  here  it  is  no  worse  than  any  other  assumption.  There¬ 
fore,  to  continue  with  our  analysis  we  will  maintain  this  assumption. 

If  Px  =  Vx  and  pt=  C2  Vt  is  assumed,  then  by  equating  and  Pt  we 
can  find  the  equivalence  between  time  and  space  variability: 
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Figure  15.  -  (Top  of  next  page).  Bias  vs.  distance  for  November  14,  1974 
at  1145  or  1215.  Three  fits  to  all  of  the  station  separation  (Table  7)  points 
are  shown  with  the  same  codes  as  in  Figure  13.  See  Table  9. 

Figure  16.  ~  (Bottom  of  next  page).  Same  as  Figure  15,  but  for  November 
15,  1974,  at  1200  or  1230.  See  Table  9. 
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P  c,  Vx 

X  =  1 


where  y  is  the  scaling  factor  (in  units  of  km/hr)  for  equating  space  and  time 
variation. 

On  November  14  y  is  (18.28/5.88)^  =  9.7km/hr  and  on  the  15th  (22. 84/3. 64)"* 

=  39.4km/hr.  Recalling  that  y  is  usually  taken  to  be  30  to  35km/hr,  we  claim 
that  the  range  of  y  found  here  on  two  consecutive  days  suggests  that  adopting  a 
single  value  for  y  or  even  a  single,  universal  weighting  system,  can  lead  to 
poorer  results  than  a  no  weighting  scheme.  Intuitively,  it  is  felt  that  y  does 
actually  vary  from  day  to  day;  there  is  no  reason  to  assume  that  the  scale  length 
of  variation  of  meteorological  parameters  over  time  and/or  space  remains  unchanged. 
In  fact,  y  is  probably  a  strong  function  of  wind  velocity  and  wind  variability. 

A  strong  steady  wind  will  undoubtedly  reduce  y  relative  to  a  no  wind  condition, 
since  it  would  take  less  time  for  a  change  in  a  meteorological  parameter  to  be 
transmitted  over  a  distance  in  a  strong  wind  than  when  it  is  calm. 

Blanco  and  Traylor  (5)  used  a  y  of  30km/hr  associated  with  a  of  0.47 
and  a  of  2.5.  While  their  y  is  roughly  similar  to  ours,  their  and 
are  an  order  of  magnitude  smaller  than  ours  on  either  the  14th  or  15th. 

Although  this  may  be  due  to  their  use  of  8-inch  howitzer  firing  tables  where 
we  used  the  155mm  tables,  it  nevertheless  emphasizes  the  danger  and  futility 
of  adopting  a  single  y  or  a  universal  weighting  system.  If  weights  are  to  be 
used,  it  is  better  to  determine  y  and  the  weights  from  available  data  each 
day. 


LEAST  SQUARES  AND  OTHERS 

A  fundamental  difference  between  our  least  squares  approach  and  the 
apparently  preferred  weighted  average  is  that  the  latter  cannot  extrapolate  to 
values  of  a  parameter  outside  the  range  already  experienced,  whereas  the 
former  can  predict  new  values.  For  instance,  suppose  that  the  measured  temp¬ 
erature  at  time  Tq-6,  Tq-4,  and  Tq-2  were,  80°,  82°,  and  84°,  respectively.  A 
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least  squares  prediction  for  the  temperature  at  Tq  would  be  86°,  but  a  weighted 
average  prediction  would  range  between  80°  and  84°  depending  on  the  weights 
employed.  Thus,  because  the  least  squares  estimator  is  not  constrained  to  a 
range,  it  can  serve  as  a  better  predictor.  This  is  a  desirable  feature, 
expecially  since  met  stations  are  often  asked  to  extrapolate  meteorological 
conditions  to  outside  the  cloud  of  stations,  to  an  artillery  battery  located 
at  or  nearer  a  battlefront. 

On  the  other  hand,  if  the  data  results  are  an  ill-conditioned  problem,  a 
least  squares  estimate  of  extrapolated  values  can  lead  to  absurd  results.  If 
two  stations  are  located  1  km  apart  and  measure  the  surface  temperature  at 
80°  and  82°,  a  least  squares  prediction  of  the  temperature  at  a  station  25  km 
away  would  be  130°.  The  weighted  average  in  this  case  would  surely  lead  to  a 
better  prediction.  However,  if  the  met  stations  are  roughly  equally  spaced, 
and  they  are  asked  to  furnish  predictions  of  parameters  not  too  far  away  in 
time  or  space  (say  not  more  than  the  average  spacing  of  the  stations  or  longer 
than  the  time  over  which  data  was  gathered),  then  we  feel  that  a  least  squares 
approach  represents  the  least  biased  method  to  adopt.  In  fact,  any  weighting 
scheme  (including  a  weighted  least  squares)  will  ultimately  compromise  the 
potential  accuracy  of  predictions  since  appropriate  weights  and  scaling  fac¬ 
tors  vary  widely  and  are  generally  unknown. 

Another  method,  the  cubic  spline  technique,  was  examined  by  D'Arcy  (lb) 
as  a  possible  interpolation/extrapolation  scheme,  but  it  was  found  "that  the 
spline  can  be  a  rather  poor  predictor"  because  in  extrapolation  "as  one  gets 
further  from  the  last  measured  point  the  slope  of  the  curve  increases  without 
bound."  Our  least  squares  approach  appears  to  offer  the  possibility  of  better 
predictions  because  unlike  the  cubic  spline,  only  one  polynomial  needs  to  be 
solved  for  each  variable,  and  not  a  polynomial  and  its  first  and  second  deriva¬ 
tives  . 

More  testing  needs  to  be  performed  to  explore  the  power  and  delineate  the 
limitations  of  the  least  squares  approach  to  extrapolating  and  interpolating 
missing  meteorological  data.  Certainly  a  safety  factor  should  be  employed  to 
prevent  absurd  values  under  unusual  conditions.  The  limits  for  extrapolation 
in  space  and  time  need  to  be  established  from  more  tesis  ur.^er  various  circum¬ 
stances.  In  addition,  the  inferences  made  concerning  space/time  variability 
should  be  confirmed  with  more  examples.  In  order  to  properly  test  a  30  km/hr 
scaling  factor,  tests  over  6  hours  of  data  should  he  compared  to  tests  over 
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(6x30=)  180  km,  whereas  we  only  tested  6  hoars  of  data  and  40  km  because  of 
the  close  station  separations.  The  choice  of  square  root,  linear,  or  power 
relations  between  bias  and  space/time  might  be  resolved  with  more  data. 
Certainly  a  greater  distance  should  be  covered  than  that  shown  in  Figures  15 
and  1 6 ,  and  a  greater  density  of  points  in  time  should  be  used  than  that  shown 
in  Figures  13  and  14. 

The  testing  that  we  have  done  on  our  least  squares  method  shows  that  it 
has  great  potential  for  extrapolating  and  interpolating  data.  At  the  very 
least  we  feel  that  the  method  is  a  viable  alternative  to  any  in  use  today. 

The  computer  code  as  it  exists  now  (to  be  documented  under  separate  cover)  is 
not  a  data  management  system  (which  is  a  necessary  and  important  part  of  the 
overall  program).  However,  it  is  flexible  enough  to  be  considered  appropriate 
to  adapt  to  a  variety  of  data  storage  systems.  Least  squares  is  an  unbiased 
and  rather  unsophisticated  approach  to  the  problem  of  missing  meteorological 
data,  but  at  the  same  time  it  is  simple  and  fairly  powerful,  easily  adaptable 
for  field  artilley  use. 


SOFTWARE  DOCUMENTATION 


1.0  INTRODUCTION 

This  documentation  pertains  to  a  least  squares  regressional  analysis  approach  to 
missing  meteorological  data,  developed  in  1981.  This  approach  relys  on  a  com* 
bination  of  horizontal  and  vertical  least  squares  fits  to  predict  missing  wind 
velocity  and  direction,  temperature,  and  pressure.  The  fitting  equations  are 
varied  by  means  of  parameters  input  from  data  cards.  Parameters  affecting  the 
fitting  equations  are  the  number  of  zones  making  up  a  met  profile,  the  number 
of  met  stations  considered,  the  number  of  different  times  considered,  and  the 
physical  location  of  the  contributing  met  stations. 

Using  data  from  seven  different  stations  collected  over  a  three  month  period, 
predictions  of  pseudo-missing  data  were  made.  Various  combinations  of  input 
stations,  times,  locations,  and  number  of  missing  layers  were  tried.  These 
predictions  were  then  used  for  a  simulated  artillery  shot  and  the  results  com¬ 
pared  with  the  actual  data  measured  at  the  point  of  interest. 

Due  to  the  fact  that  this  approach  uses  all  available  data  as  input,  only 
missing  data  from  one  site  at  a  time  may  be  predicted.  Obviously,  well 
behaved  input  will  provide  more  accurate  results,  which  will  enhance  the 
value  of  the  predictions  as  input  to  later  predictors.  Hence,  it  is  well 
worth  the  effort,  if  possible,  to  analyze  and  discard  rough  or  inaccurate 
inputs . 

At  this  stage  of  development,  this  software  would  be  a  valuable  "front-end" 
tool  for  artillery  meteorological  units. 


2 . 0  INPUT 


Card  Input 

Cols . 

CARD  ONE: 

Format 

Name 

Description 

1-3 

13 

NS 

Number  of  sites  input 

4-6 

13 

NT 

Number  of  times  input 

7-9 

13 

NZ 

Number  of  zones 

CARD  TWO: 

1-8 

F8.0 

CX 

X  coordinate  of  sta. 

9-16 

F8.0 

CY 

Y  coordinate  of  sta. 

17-24 

F8.0 

CZ 

Z  coordinate  of  sta. 

25-32 

F8.0 

THR 

Hour  of  met  message 

33-40 

F8.0 

TMIN 

Minute  of  met  message 

CARD  THREE: 

1-12 

2(A6) 

SITE 

Two  word  array  containing 
site  ID,  time  and  date  to  be 
read  from  mass  storage  device 

CARD  FOUR: 

1-2 

12 

NZONS 

Number  of  layers  in  profile 
defined  by  card  three 

CARD  FIVE: 

1-3 

13 

NINPRO 

Number  of  layers  in  each  met 
profile 

4-6 

13 

MISLYR 

Number  of  missing  layers  to 
be  predicted 

Card  two  is  repeated  for  each  profile  input.  The  values  on  the  first  card  two 
become  the  values  of  the  origin.  The  values  on  the  last  card  two  pertain  to  the 
station  with  missing  data. 
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In  che  present  configuration,  cards  three  and  four  are  repeated  as  a  pair.  One 
pair  for  each  card  two.  This  will  change  as  the  software  is  implemented  with 
different  data  input  devices  and  formats. 

Cards  one  and  five  are  both  read  in  from  the  main  program  and  neither  their 
number  nor  position  should  change. 

Currently,  test  data  is  read  from  mass  storage  using  logical  unit  seven.  This 
data  is  in  a  Fortran  formatted  file.  Each  profile  consists  of  twelve  lines 
of  data,  a  header  line  in  the  same  format  as  input  card  three,  ten  data  lines 
consisting  of  layer  number,  wind  direction,  wind  speed,  temperature,  and 
pressure.  The  data  lines  are  formatted,  12,  13,  13,  15,  14,  and  the  last 
data  line  is  followed  by  a  line  containing  99,  which  denotes  the  end  of  profile. 


3.0  OUTPUT 

Program  output  is  currently  directed  to  a  line  printer,  logical  unit  six. 
Output  consists  of  two  parts,  the  first  part  is  merely  an  echo  of  the  input 
coordinates  and  times  printed  in  a  5(F8.0)  format.  The  second  part  of  the 
output  is  the  met  profile  of  interest,  with  the  actual  data  as  far  as  it  were 
available  and  predictions  finishing  out  the  profile.  These  data  are  printed 
in  the  following  format: 

Layer  Number  (14),  Wind  Direction  (F10.3),  Wind  Speed  (F10.3), 
Temperature  (F10.1),  and  Pressure  (F10.0). 
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4.0  OPERATING  INSTRUCTIONS 

The  program  will  be  provided  on  two  media,  punched  card  and  nine~track 
magnetic  tape. 

Punched  card  decks  contain  all  necessary  control  cards  to  compile  and  assemble 
the  routines  into  an  absolute  element.  Identical  source  language  and  control 
cards  are  on  tape  in  a  file  called  8102*PREDICT. . 

Source  programs  are  in  FORTRAN  V.  Only  one  change  is  necessary  to  convert  to 
ASCII  Fortran.  The  four  @F0R,IS  cards  need  to  be  changed  to  @FTIN,IS.  All  other 
cards  and  procedures  remain  the  same. 
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PROGRAM  VARIABLES 


Variables  In  Common 


CX,  CY,  CZ  » 
INDEX  » 

JZONS  ■ 

NA,  NB,  NQ,  N,  M 
NINPRO  « 

NSITES  * 

NZONS 

Range  &  Cross  * 
SE 

TEMP  &  PRES 


THEATA 

THR  &  MIN 

TRANGE 

TCROSS 

TPRES 

TTEMP 

VEL 


X,Y,Z  coordinates  of  each  net  station 

An  error  indicator  from  LSTSQR  not  currently  being 

output 

The  total  number  of  layers  input 
*  Input  arguments  to  LSTSQR 
The  number  of  layers  in  the  net  profile 
The  number  of  profiles  input 

Input  value  giving  the  number  of  layers  in  the  input 
profile 

Measured  X  and  Y  displacements  at  each  level 
Output  standard  error  of  single  equation  from 
subroutine  LSTSQR 

Measured  temperature  and  pressure  values  at  each 
level 

The  input  wind  direction  in  miles 

The  hour  and  minute  time  log  for  each  profile 

Predicted  value  of  range 

Predicted  value  of  cross 

Predicted  value  of  pressure 

Predicted  value  of  temperature 

Input  wind  velocity  in  knots 

Cummulative  X,  Y,  Z  and  time  balloon  displacement 
computed  from  the  origin 


Local  Variables  In  Main 


A 

CC,  OMCC,  CE 

DELT 

DELZ 

LSTKWN 

MISLYR 

NS 

NT 

NZ 

RQ,  OMC,  C 


Input  array  of  coefficients  for  LSTSQR 

Intermediate  arrays  used  to  compute  predictions 

Difference  between  T  (Istkwn)  and  T  flstkwn-1) 

Difference  between  Z  (Istkwn)  ar.d  Z  ( Istkwn*  1) 

INDEX  at  last  known  X,  Y,  Z  and  time 

Number  of  missing  byers 

Number  of  different  sites  input 

Number  of  different  times  input 

Number  of  zones  input 

Output  arrays  from  LSTSQR 


III.  Local  Variables  in  LSTSQR 


A  »  Coefficients  of  observational  equations 

C  *  Computed  value  for  all  observatioal  equations 

#MC  *  Observed  -  Computed  for  all  observational  equations 

Q(N,J)  *  X(J)  the  unknowns 

Q ( I , J )  *  The  weights  of  the  unknowns 

Q(N,N)  *  The  sum  of  the  (0-C)  squared 

R(I,J)  *  The  coefficients  of  the  normal  equations 

R(N,J)  *  The  standard  error  of  the  unknowns 

R(N,N)  *  The  sum  of  the  (f-C)  squared,  calculated  from  0MC  and 
used  to  calculate  SE 

SE  *  The  standard  error  of  a  single  equation  of  unit  weight 


IV.  Local  Variables  in  READER 

Header  *  Two  word  array  with  site  id,  date  and  time,  read 
from  mass  storage 

SITE  *  Two  word  array  containing  site  id,  date  and  time  of 
interest 

ITHETA  *  Temporary  variable  used  to  read  integer  wind  direction, 
units  are  10's  of  mils 

IVEL  =  Temporary  variable  used  to  read  integer  wind  velocity  in 
whole  knots 

ITEMP  »  Temporary  variable  to  read  integer  temperature,  units  are 
tenth's  of  degrees 

IPRES  *  Temporary  variable  used  to  read  iiteger  pressure  in  millibars 


V.  Local  Variables  in  C0EF 

All  variables  except  indices  in  C0EF  are  contained  in  common. 
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nnonnnoonnnnnnononnnoonno 


@F0«,1S  main 

COMMON  / A B /  $E,CX'13) ,CY(13) ,C2( 13) ,THR{13)  ,TMIN(  13) , 

1  THETA (  340  ) ,VEL(340) .NINPRO, 

1X(340),Y(340),Z(340),T(3401, 

IRANGEf  340) .CROSS (  340  ) .TEMP (  340  )  .PRES 040), 

1TRANGE(26) ,TCR0SS(26)  , 

1TTEMP ( 26 ) ,TPRES(26) , 

1 N Z 0 N S ( 13)  .JZONS.NSITES 
1, NA.NB, NO, N.M, INDEX 

DEFINITION  OF  VARIABLES  IN  COMMON 

SE-OUTPUT  STANDARD  ERROR  OF  SINGLE  EQUATION 
FROM  SUBROUTINE  LSTSQR 

CX.CY.AND  CZ*X,Y,Z  COORDINATES  OF  EACH  MET  STATION 
THR  AND  TM I N *  THE  HOUR  AND  MINUTE  TIME  TAG  FOR  EACH  PROFILE 
THETA-  THE  INPUT  WIND  DIRECTION  IN  MILLS 
VEL-  THE  INPUT  W’ND  VELOCITY  IN  KNOTS 
NINPRO-  THE  NUMBER  OF  LAYERS  IN  THE  MET  PROFILE 
X.Y.Z.T-  CUMULATIVE  X.Y.Z.AND  TIME  BALLOON  DISPLACEMENTS 
COMPUTED  FROM  THE  ORIGIN 
RANGE  AND  CROSS-  KNOWN  X  AND  Y  DISPLACEMENTS  AT 
EACH  LEVEL 

TEMP  AND  PRES-  KNOWN  TEMPERATURE  AND  PRESSURE  VALUES 
AT  EACH  LEVEL 

TRANGE  .TCROSS  TTEMP  .AND  TPRES*  PREDICTED  VALUES  OF 
RANGE, CROSS, TEMP, AND  PRES 
NZONS*  INPUT  VALUE  GIVING  THE  NUMBER  OF  LAYERS  IN 
THE  INPUT  PROFILE 

JZONS*  THE  TOTAL  NUMBER  OF  LAYERS  INPUT 
NSITES-  THE  NUMBER  OF  PROFILES  INPUT 
NA,NB,NQ,N,M-  INPUT  ARGUMENTS  FOR  LSTSOR 
INDEX-  AN  ERROR  INDICATOR  FROM  LSTSQR  WHICH  IS  NOT 
CURRENTLY  BEING  OUTPUT 

DOUBLE  PRECISION  A ( 340 , 33 ) , R (  33 , 33  ) . Q (  3 3 , 3 3  ■  .CMC ( 340 ) , C ( 340 ) 
DOUBLE  PRECISION  CG 

DOUBLE  PRECISION  CC ( 340 ) ,OMCC < 340 ) , CE ( 340 ) 

C 

C  LOCAL  VARIABLES 
C 

C  A-  INPUT  ARRAY  OF  COEFFICIENTS  FOR  LSTSQR 
C  R  ,Q  ,OMC .AND  C  -  OUTPUTS  FROM  LSTSQR 

C  CC.OMCC.CE-  , INTERMEDIATE  ARRAYS  USED  TO  COMPUTE  PREDICTIONS 
C  NS-  NUMBER  OF  DIFFERENT  SITES  INPUT 
C  NT-  NUM8ER  OF  DIFFERENT  TIMES  INPUT 
C  NZ-  NUMBER  OF  ZONES 
C  MISLYR-  NUMBER  OF  MISSING  LAYERS 
C  LSTKWN-  INDEX  OF  LAST  KNOWN  X.Y.Z.AND  TIME 
C  CELT-  DIFFERENCE  BETWEEN  T (LSTKWN )  AND  T(LSTKWN-l) 

C  OELZ-  DIFFERENCE  BETWEEN  Z ( LSTKWN )AND  Z(LSTKWN-l) 

READ ( 5 , 98  )NS  , NT  ,  NZ 
98  FORMAT (3(13)) 

NA- 340 
INDEX-1 
NB-33 
NQ-33 

NSITES«NS*NT 
CALL  READER 
CALL  COEF 

REA0(5,199)NINPR0, MISLYR 
199  FORMAT (13,13) 

KNOWN- ( NS  I TES- 1  )*NINPRO 
ICUT-NINPRO-MIELYR 
1 1  - 1  CUT 
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li>lKWN--KNOWN+l<  It T 

NN  *0 
NM  *  0 
l  I  = I  I +K 

NN*( ( I0PT-1 )*340) 

NZZ  *N  Z 
NSS*NS 
NTT*NT 

00  600  I  OPT*  1 , 3 
NZ=NZ2 
NS*NSS 
NT=NTT 

I F  f I  OPT . GE . 3  )G0  TO  2000 
C  BEGIN  HORIZONTAL  FIT 
1*1 
0*1 

M*NS*NT - 1 
GO  TO  2 

3  J*J  +  1 
GO  TO  2 

4  1=1+1 

2  N=1+2*(NS-I  )  +  ( NT- J  ) 

I F ( M  . GT . N  ) GO  TO  1000 
IF( (NS.EQ.l ) .AMD. CNT.EQ.2)  )G0  TO  1000 
IF(N.LE.l)Gt)  TO  2000 
I F ( (NT-J) .GT. (NS-I  )  )G0  TO  3 
GO  TO  4 
1000  NZ*0 

GO  TO  2200 

C  BEGIN  VERTICAL  FIT 
2000  I F ( NZ . GT .  5  )  NZ-5 

IF{  I  OPT. EO. 4)  GO  TO  2010 
GOTO  2080 
2010  NS  *  1 
NT»0 
N  Z  *  1 

2080  1-1 
J  =  0 
K  *0 

M=NZZ*NS*NT-MISLYR 
3100  J=J  +  1 

3200  N=1+2*(NS-J)+NT-J+NZ-K 

I F ( M . GT .  N  )G0  TO  2100 
IFfN.LE.  (NZ-K+1 ) ) GO  TO  4000 
IF( (NT-J)  .GT.  (NS-I  ) )  GO  TO  3100 
1*1  +  1 

GO  TO  3200 
4000  K*K+1 

GO  TO  3200 
2100  NZ«NZ-K 
2200  NT*NT- J 
NS-NS-I 
N*N  +  1 

IF ( IOPT  .EQ .  3 )G0  TO  700 
GO  TO  500 

C  500  COMPUTES  THE  A  ARRAY  FOR  HORIZONTAL  FIT 
500  DO  117  LUP-l.MISLYR 
1 1 *1 1 +LUP 
NM-0 

NN*( ( I  OPT  - 1  ) *340 ) 

C 

C  OECISION  FOR  HORIZONTAL  OR  VERTICAL  FIT 
C 

IF(I0PT.LE.2)L0P$IZ«N$ITES 

IF(IOPT.E0.3)LOPSIZ-JZONS 

LST-LSTKWN-1 
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DLL  T  *  1  'Li'!  KWN  1  -  T  t  L  S  T  > 

OELZ*ZaSTKWN)-ZfLST' 

MM*NN 

00  5  1  *1 , LOPS  r  z 
J *  1 1 *NN 
J J* I I+NM 

IF(I.GT.«)X(JJ)-X(LSTKWM) 

IF ( I .GT.M)YfOJ)-Y(LSTKWN) 

I  F  (  I  ,GT?M)T(  JJ)«T(LSTKWN)*(DELT*LUP) 

I F  ( I  .GT.M)Z(  JJ)-ZaSTKWN)*(OELZ*LUP) 

A (  I  ,1)-1. 

IF(NS.EQ.O)GO  TO  11 

00  10  K-l.NS 

KK-K+l 

KKK*KK+NS 

A ( 1 ,  KK  )  *X ( J J  )  **K 

10  A ( I , KKK ) « Y ( J J )  **K 

11  CONTINUE 

IF(NZ. £0.0)60  TO  21 
DO  20  K«1,NZ 
N Y *K  + 1  +  ( 2*N S  ) 

20  A( I , N  Y ) *  Z ( J J ) **K 

21  CONTINUE 

I F ( NT . EO . 0 )G0  TO  31 
DO  30  K  *  1 ,  NT 
I  N*K  + 1  ■*■  (  2*NS  )  *N  Z 

30  A( I  ,1N)«T( JJ)**K 

31  CONTINUE 
INO*NN+I 

A (  I  ,  N  )  ’RANGE ( J  ) 

NN  *NN  +  N  I  NPRO 
NM=NM+N I NPRO 
5  CONTINUE 

C 

C  CALL  LSTSQR  TO  PRF.OICT  MISSING  WIND  LAYER 
C 

CALL  LSTSQR(A,r.,q,OMC,C) 

CG*0 
INOX  *N- 1 

DO  188  I T » 1 , 1 NDX 
CG*CG+Q(N, IT) *A (NSITES.IT) 

188  CONTINUE 

I F ( MM  . GT . 0  )G0  TO  113 
TR ANGE ( 1 1  )  =CG 
GO  TO  117 
113  CONTINUE 

TCROSS ( 1 1  ) -CG 

117  1 1  *  I  CUT 
111  CONTINUE 

DO  118  I - 1 , 1  CUT 
K  *  I +KNOWN 

TRANGE { I ) -RANGE ( K  ) 

118  TCROSS ( I ) "CROSS ( K  ) 

GO  TO  12 

121  CONTINUE 

700  NN*  ( ( I  OPT- 1  ) *340 1 

IF(I0PT.LE.2)t.0PSIZ»NSITES 

IF(10PT.EQ.3)L0PSIZ»JZ0NS 

C 

C  COMPUTES  A  ARRAY  FOR  VERTICAL  TEMPERATURE  FIT 
C 

DO  15  I  *  1 .LOPS I Z 
I F ( I .GT.M)X(I ) -X (M ) 

IF ( I . GT .M  )  Y ( I  )  * Y ( M ) 
IF(I.GT.H)T(I>T(I-1)*0ELT 
1F(I.GT.M)Z(I  WM-n+DElZ 
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A  (  I  ,  1 )  ■  1 . 

I F ( NS . EQ . 0 'GO  TO  all 

DO  410  K«1,NS 

KK«K+1 

KKK-KK+NS 

A ( I ,KK )»X( I  }**K 

410  A( I ,KKK)»Y( I )**K 

411  CONTINUE 
IF(NZ.EQ.O)GO  TO  421 
00  420  K*1,NZ 
NY*K+1+ ( 2*NS ) 

420  A ( I , N Y ) »Z ( I ) **K 

421  CONTINUE 

I F ( NT . EQ . 0 )G0  TO  431 
DO  430  K«1,NT 
I N=K  +  l  +  ( 2*NS )  *NZ 

430  A ( I  ,  I  N  )  =T (  I  )  **K 

431  CONTINUE 
I ND=NN+ I 

A ( I , N ) *RANGE ( I  NO  ) 

15  CONTINUE 
C 

C  CALL  LSTSQR  TO  PREDICT  MISSING  TEMP  LAYERS 
C 

CALL  LSTSQR(A,R,Q,OMC,C) 

L*1 

MM=M+1 

IA=KN0WN+NINPR0 
DO  88  I  =MM  ,  I A 
INOX=N- 1 

00  187  IT-1, INOX 
CC(L)-CC(L)+Q(N,IT)*A(I,IT) 

187  CONTINUE 

OMCC(L)=A( I  ,N)-CC(L) 

L  =  L  +  1 

88  CONTINUE 

NM= JZONS-M 
MN=KN0WN+1 
MM=MN 

DO  33  1=1, 1  CUT 
K= I +KN0WN 
TTEMP ( I ) =TEMP ( K  ) 

33  CONTINUE 

L*1 

IB=NINPR0-MISLYR+1 
00  32  I  *  I B , N I NPRO 
TTEMP ( I  )=CC ( L  ) 

L  *L  +  1 

32  CONTINUE 

12  CONTINUE 

600  CONTINUE 

M«KNOWN+(NINPRO-MISLYR) 

N«3 

C 

C  COMPUTES  A  ARRAY  FOR  VERTICAL  PRESSURE  FIT 
C 

DO  211  1  =  1,  JZONS 
A  ( I  ,  1 )  « 1 .  DO 
A{  I  ,  2  )  =  Z ( I  ) 

IF ( PRES ( I ) . GT .0 . ) A ( 1 , 3 ) -OLOG ( PRES n ) ) 

I F ( PRES ( I ) . EO . 0 . ) A( 1 , 3 ) =0 . 

211  CONTINUE 
C 

C  CALL  LSTSQR  TO  PREDICT  MISSING  PRES  LAYERS 


i  &  J 


CALL  LSTSQR(A,R,n,UMC  ,C  ) 

MM  *M+ 1 

DO  222  I *MM , J ZON S 

222  C(n*Q(N,l)*AU,l)*0(N,2!*A(l,2' 

DO  987  I ■  1 , JZONS 
CE ( I ) *DEXP ( C (  I  )  ) 

987  CONTINUE 

MN*KN0WN+1 

MM*MN 

DO  338  1*1  ,  1  CUT 
TPRES ( I  ) -PRES (MN  ) 

MN*MN+1 

338  CONTINUE 
IB*NINPRO-MISLYR+l 
DO  339  I  *  I B , N I NPRO 
TPRES ( I ) *CE (MN  ) 

MN-MN+1 

339  CONTINUE 
C 

C  CONVERT  PRESDICTED  VALUES  BACK  TO  STANDARD  MET  UNITS 
C 

DO  341  I = 1 , N I NPRO 
TDIR=ATAN(TRANGE( I ) /TCROSS (  I ) ) 

TD IR=( TDIR/ (2*3. 14159 ) )*6400. 

TDIR=4800-TDIR 

IF (TDIR .LT.O. )TDIR=TDIR+6400. 

TVEL  =  (SQRT( (TRANGE{ I ) **2 )  +  ( TCROSS ( I ) **2 )))/(. 5 1444444*60  .  ) 
WRIT£(6, 1357)1, TDIR, TVEL,TTEMP(I),TPRES{I) 

1357  FORMAT ( IX, 14 tF10.3,F10.3,F10.1,Fl 0.0) 

341  CONTINUE 
STOP 
END 
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OFOR.IS  LSTSOR 

SUBROUTINE  L STSQR ( A , R , 0  , OMC  , C  I 

COMMON  /AS/  SE,CX(13),CY(13),CZ(13).THRn3),TMIN{13), 
1THETA(340),VCL(340),NINPR0, 

1X(  340) ,Y ( 340) , Z ( 340) ,T (  3401  , 

IRANGEf 340), CROSS (340), TEMP (340), PRES (340), 

1TRANGE ( 26 ) ,TCR0SS(26), 

1TTEMP(26),TPRES(Z6), 

INZ0NS( 13  )  ,J20NS,NSITES 
1 ,NA,NB,NQ,N,M, INDEX 

DOUBLE  PRECISION  A ( 340 , 3 3 ) , R ( 3 3 , 3 3 ) , Q ( 3 3 , 3 3 ) , OMC < 340 ) , C ( 340 ) 

C  LEAST  SQUARES  SOLUTION 

C  INPUT  ARGUMENTS  ARE  A , NA , NQ , N , M . 

C  OUTPUT  ARGUMENTS  ARE  R , 0 , SE , OMC , C  . 

C 

C 

C  INDEX  IS  THE  ERROR  COMPUTATIONAL  SWITCH 
C  IN0EX=0  MEANS  SUCCESFUL  TERMINATION 
C  I NDEX  =  2  MEANS  BAD  N,  OR  M 

C  I NDEX  =  3  MEANS  Q(I,I)IS  LESS  THAN  OR  EQUAL  TO  ZERO 

C  TO  SUPPRESS  THE  ERROR  MESSAGE  I F ( Q ( I  ,  I )  I S  LESS  THAN  0..  ENTER  LSTSORS 
C  WITH  INDEX  =0.  ,0R  1  . 

C  I F  Q (  I  ,  1  )  L  E  0  THE  COEFFICIENTS  OF  THE  NORMAL  EQUATIONS  WILL  STILL 
C  HAVE  BEEN  CORRECTLY  FORMED  AND  RETURNED.  LIKEWISE  FOR  ALL 
C  PREVIOUS  0(1,1)  S  AND  THEIR  0(1,  J)S  WHERE  J  GT  I. 

C 

C 

C  NOTATION 
C 

DUM  =  1 

N=NUMBER  OF  UNKNOWNS  °LUS  1 

M -NUMBER  OF  OBSERVATIONAL  EQUATIONS  OF  CONDITION, 

A ( K  ,  J  )  =  COEFFICIENTS  OF  THE  OBSERVATIONAL  EQUATIONS  OF  CONDITION 

A  (  K  ,  1  )  *X  (  1  )  +  A  (  K  ,  2  )  *X  (  2  )  ♦ . +  A  (  K  ,  N  -  1  '  *  X  ■'  N  -  1  ' 

=  A  (  K  ,  N  ) 

WHERE  K  =  1 ,  2  ,  3 . M 

R (  I  ,  J  )  =  COEFFICIENTS  OF  THE  NORMAL  EOUATIONS,  WHERE 

1  =  1,2,  ...  ,N-1  AND  J=I , 1+1 , 1+2 . N, 

R ( N  ,  J  )  =  STANDARD  ERROR  OF  THE  UNKNOWNS  ,  J*  1  ,  2 . N-l, 

0 ( N  ,  J  )  =  X(J),  THE  UNKNOWNS  J  =  1  ,  2 . N-l 

Q(I, J)=  THE  WEIGHTS  OF  THE  UNKNOWNS  J  =  1  ,  2 . N-l, 

AND  1=J  ,  J+l ,  .  .  .  ,N-1 

OMC ( K ) =  OBSERVEO-COMPUTED  FOR  ALL  03SER VAT  I ONAL  EQUATIONS, 

C ( K ) =COMPUTED  VALUES  FOR  ALL  03SE R V  AT  I  ON AL  EOUATIONS  OF  CONDITION 
K-1,2,3 . M 

SE=STANDARD  ERROR  OF  A  SINGLE  EQUATION  OF  UNIT  WEIGHT, 

=  THE  SQUARE  ROOT  OF  THE 

=  SUM  OF  ALL  ( 0-C )  **2  DIVIDED  BY  (M-fN-H) 

(SE  MULTIPLIED  BY  THE  SQUARE  ROOT  OF  THE  SUM  OF  THE 
SQUARES  OF  THE  WEIGHTS  OF  AN  UNKNOWN  IS  THE  UNKNOWN  S 
STANDARD  ERROR  ) 

R ( N  ,  N  )  *  SUM  OF  (0-C  )**2  CALCULATED  FROM  0MC(K),K«1,2 . M 

(AND  USED  TO  CALCULATE  SE )  , 

Q ( N , N )  =  SUM  OF  THE  (0-C)**2  CALCULATED  FROM  THE  IDENTITY  -- 

( 0-C ) **2  =SUM  OFA(K,N)*A(K,N)  WHERE  K-1,2 . M 

MINUS  THE  SUM  OFQ ( K  , N  ) *Q ( K  , N  )  WHERE  K«1 ,2  ,  .  .  .  ,N-  1 
(NOT  USED  TO  CALCULATE  SE ) 

NA  =THE  NUMBER  OF  ROWS  IN  THE  A  ARRAY  AS  DIMENSIONED  IN  THE 
MAIN  PROGRAM  (M  MUST  BE  LE  TO  NA ) 

NQ*  THE  NUMBER  OF  ROWS  IN  THE  N  AND  Q  ARRAYS  AS  DIMENSIONED 
IN  THE  MAIN  PROGRAM (  THE  R  AND  Q  ARRAYS  SHOULO  BE  OF  THE 
SAME  SIZE  AND  N  MUST  BE  LE  TO  NO) 
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c 

NM 1 =N -  1 

IF((M.GE.NM1).OR.(M.LE.NA).OR.(N.LE.NQ).OR.(N.GE.2))GO  TO  20 
WRITE (6,2) 

2  FORMAT ( 10X , ' NO  OF  OBSERVATIONS  LT  NO  UNKNOWNS  OR  TO  HAVE  EXCEED', 

1 ' ED  DIMENSIONS  OR  NO  OF  UNKNOWNS  LT  1') 

I  NOE  X  =  2 
RETURN 
C 

C  CALC.  COEFF.  OF  NORMAL  EQNS  .  ,  R. 

20  DO  30  1  =  1, N 

I K  =  I 

DO  30  J= I K  ,  N 
R ( I , J ) »0 . DO 
DO  30  K  =  1 ,  M 

30  R(I,J)-R(I.J)+(A(K,I)*A(K,J)) 

C 

C  CALCULATE  TRIANGULAR  SQUARE  ROOT  OF  R  CRACOVIAN  AND  THE  RECRIP.  FO 
C  ITS  DIAGONAL  ELEMENTS  (EXCEPT  IF  I*N  DONT  TAKE  ITS  RECRIP. 1 
C 
C 

C  EVALUATE  0(1,1). 

DO  120  1  =  1, N 
Q ( I  ,  I  ) «R( I  ,  I ) 

I M 1  =  I  - 1 

IF ( I M 1  .EQ.O)GO  TO  41 
DO  40  K  =  1 , 1 M 1 

40  Q(I,I)=Q(I,I)-(Q(K,I)*Q(K,I)) 

41  CONTINUE 

I F ( I  .  EQ .  N  )G0  TO  130 

C 

C  ERROR  CHECK  POINT  2 

C 

C 

I F ( Q ( I ,  I  )  .  GT  .  0  .  DO  )G0  TO  80 
50  IF(  INDEX. LT.2)G0  TO  75 

60  WRITE (6,61) 

61  FORMAT ( 10X  , ' NEG  ARG  IN  SORT') 

WRITE (6, 70)I,I,Q(I,I) 

70  F0RMAT(1X,/,10X, 'THE  SQUARE  OF  THE  ',12,'  ’,12,'  ELEMENT  OF  THE', 

1 ' TR I  ANGULAR  SQUARE  ROOT  OF  THE ' , / , ' 1  OX . MATR I X-CRACOV I  AN  CONTAIN', 
2'ING  THE  COEFFICIENTS  OF  THE  NORMAL  EQUATIONS  IS  '.E15.6) 

75  INDEX=3 

RETURN 
C 
C 

80  Q(I,I)=1.D0/DSQRT(Q(I,I)) 

90  IP1=I+1 

IF ( I P 1 . GT . N )G0  TO  130 
C  EVALUATE  Q(I,J)FOR  ALL  J  GT  I 
DO  120  J  *  I P 1 , N 
Q(  1 , J ) =R ( I , J  ) 

I F ( IM1 . EQ . 0 )G0  TO  101 
DO  100  K*1 ,  IM1 

100  Q(I,J)*Q(I,0)-(Q(K,I)*Q(K,J)) 

101  CONTINUE 

120  Q( I , J)*Q( I , J)*Q( 1,1) 

C 

C 

C  EVALUATE  Q(I,J)  FOR  ALL  J  LT  I. 

130  DO  160  J* 1 , N 

JP1-0+1 

IF(JPl.GT.N)  GO  TO  160 
DO  150  I  * JP1 , N 
Q  ( I  ,  J  )  *0 . DO 

TM1  «  T  .  1 
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I F ( IM 1 . EQ . 0 )G0  TO  141 
DO  140  K=J, IM1 

Q(I,J)'Q(I,J)+(Q(K,I)*0(K,J)) 

140  CONTINUE 

141  CONTINUE 

IF( I  .EQ.N  ) GO  TO  160 
QU,j)=-(Q(i,j))*Q(i,n 
150  CONTINUE 
160  CONTINUE 

C  CALC  COMPUTED  VALUES ,( 0-C ) .  AND  STD  ERROR  OF  EACH  EQN  OF  UNIT  WT. 

R  (  N  ,  N ) =0 . DO 
C 

DO  180  K  =  1 , M 
C  ( K  )  =0 . 

DO  170  J *  1 , NM 1 

170  C(K)=C(K)+(A(K,J)*Q(N,J)  ) 

OMC(K)=A(K,N)-C(K) 

180  R  (  N  ,  N  )  *R  (  N  ,  N  )  +  ( OMC ( K  )  *OMC ( K  )  ) 

C 

SE  =  DS0RT(R(N,N)/(M-NM1)) 

NM2  =  NM 1 - 1 

C  CALCULATE  STD  ERROR  OF  UNKNOWNS 
DO  200  J  =  1 ,  NM2 
R ( N , J ) =0 . DO 
DO  190  I = J , NM1 

190  R(N,J)*R(N,J)  +  (Q( I ,J)*Q{  I  ,J)  ) 

200  R(N,J)=SE*DSQRT(R(N,J)  ) 

C 

R(N,NM1)=SE*Q(NM1,NM1) 

C 

I NDEX=0 

RETURN 

END 
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0FOR.IS  REAO 

SUBROUTINE  REAOER 

COMMON  / A  8 /  SE,CX( 13) ,C Y ( 1 3 1 , CZ ( 1 3 ) , THR ( 13  )  ,TMI  N  f 1 3 )  , 

1  THETA ( 340 ) ,VEL(  340} .NINPRO, 

1 X (  340  ) , Y (  340 ) ,  Z  (  340) , T ( 340 )  , 

1RANGE(  340) .CROSS (  340) , TEMP (  340) , PRES (340) , 

1TRANGE(26) ,TCR0SS(26)  , 

1TTEMP(26).TPRES(26). 

1NZ0NS(13) .JZONS.NSITES 
1 ,NA,NB,NQ,N,M, I  NOE  X 
DIMENSION  HEADER ( 2 )  ,  S I TE ( 2  ) 

C 

C  LOCA?  VARIABLES 
C 

C  HEAOER-  ARRAY  WITH  SITE  ID, DATE  AND  TIME  READ  FROM 
C  _MASS  STORAGE 

C  SITE  =  ARRAY  CONTAINING  SITE  ID, DATA, AND  TIME  OF  INTEREST 
C  ITHETA, IVEL , ITEMP.AND  IPREA  ARE  TEMPORARY  VARIABLES 
C  ALLOWING  THE  READING  OF  INTEGER  VALUES 
C 

JZONS-1 

DO  11  L-l.NSITES 

REAO( 5,201 )CX(L) ,CY(L) ,CZ(L) ,THR(L)  ,TMIN(L) 
WRITE(6,201)CX(L),CY(L) ,CZ(L),THR(L),TMIN(L) 

C 

C  CONVERT  FEET  TO  METERS 
C 

CX(L)-CX(L)*.3048 
CY ( L ) -C Y ( L  )  *  . 3048 
CZ ( L ) =CZ ( L ) * . 3048 
11  CONTINUE 
201  FORMAT ( 5 ( F 8 . 0 )  ) 

200  FORMAT (12) 

DO  2  IK-1, NSITES 
REWIND  7 

READ(5,501)SITE(1),SITE{2) 

501  F ORMAT ( 2 ( A6  )  ) 

C 

C  SEARCH  MASS  STORAGE  FOR  SITE 
C 

1  READ(7,501)HEADER(1),HEADER(2) 

IF( (HEADER (l).NE.SITE(l) ) . OR . ( HEADER 1 2 ) . NE . S I TE ( 2 ) )  ) 

1G0  TO  1 

READ(5,200)NZ0NS(IK) 

NK=NZONS( IK) 

C 

C  READ  LAYERS 
C 

DO  3  I«1,NK 
IF(NK.EQ.O)GO  TO  3 

RE AD (7, 300) I  LEV, ITHETA, IVEL, ITEMP,  I  PRES 
THETA( JZONS  ) « ITHETA 
VEL( JZONS ) - I VEL 
TEMP ( JZONS  )  ■ I  TEMP 
PRES( JZONS)«IPRES 
300  F0RMAT(I2,2(I3),I5,I4) 

JZONS- JZONS+1 
3  CONTINUE 

2  CONTINUE 
JZONS-JZONS-1 
RETURN 
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@FOR,IS  COEF 


SUBROUTINE  COEF 

COMMON  / AB /  SE,CX(13),CY(13),CZfl3),THR(13),TMIN(13), 
1TH£TA(340),VEL(340),NI6PR0, 

IX(340),r(340),Z(340),Tf340), 

1RANGEI340) .CROSS (340) .TEMP (340), PRES (340) , 

1TRANGE(26) ,TCROSS( 26)  , 

1TTEMP ( 26  )  ,  TPRES ( 26  )  . 

IN  ZONS (13) .JZONS.NSITES 
1 ,NA,NB,NQ,N,M, INDEX 
C 
C 

C  ALL?  VARIABLES  EXCEPT  INDICES  ARE  IN  COMMON 
C 

C  CHNGE  VELOCITY  TO  METERS/MTNUTE  AND 
C  UNSCALE  TEMPERATURE  ANO  PRESSURE 
C 

DO  10  1=1, J  ZONS 

VEL ( I )=VEL( I )*( .51444444*60.  ) 

TEMP ( I ) =TEMP ( I ) / 1 0 . 

10  THETA( I )=(THETA( I ) )*10. 

C 

C  CONVERT  FROM  NAVIGATIONAL  COORDINATE  SYSTEM  TO 
C  MATHMATICAL  COORDINATE  SYSTEM 
C 
C 

I  Z  =  1 

DO  20  K=1 .NSITES 
NZ  =  NZONS ( K  ) 

DO  20  I K  =  1  ,  N  Z 
I  =  1 Z 

IF(THETA( I ) .GT.4800.  )  GO  TO  30 
ANGL=4800. -THETA( I  ) 

GO  TO  40 

30  ANGL=(4800.-THETA(I))+6400 

40  CONTINUE 

ANGL=(ANGL/6400. ) * ( 2 . *3 . 141 59 ) 

C 

C  RANGE  IS  THE  NORTH, SOUTH  COMPONENT  OF  THE  WIND 
C  NORTH  BEING  POSITIVE 
C 
C 

C  CROSS  IS  THE  EAT, WEST  COMPONENT  OF  THE  WIND 
C  EAST  BEING  POSITIVE 
C 

RANGEU  )  =  SIN(ANGL)*VEL(I) 

CROSS ( I )=COS(ANGL)*VEL( I ) 

IF(  IK.EQ.DGO  TO  21 
I F ( IK  . GT . 4  )G0  TO  33 
C 

C  X  AND  Y  ARE  DISPLACEMENTS  COMPUTED  USING  VELOCITY  ANO  TIME 
C  Z  AND  T  ARE  CUMULATIVE  HEIGHT  AND  TIME  RESPECTIVELY 
C  Z  IS  COMPUTED  ASSUMING  A  CONSTANT  RISE  RATE  c  300  METERS  PER  SECOND 
C 

IF(IK.EQ.2)X(I)=CR0SS(I)*. 33333 
I F ( IK . EQ . 2 ) Y ( I )«RANGE( I  )*. 33333 
IF(IK.EQ.2)Z(I)«100. 

IF(IK.EQ.2)T(I)«. 33333 

IF(IK.EQ.3)X(I>«(CROSS(I-l)*.33333)*(CROSSU  )*.5) 

I F  (  IK  .  EQ .  3 )  Y  ( I  )«(  RANGE  {  I  -  P*.  33333  )*  (RANGE  f  I  '*.5) 

IF ( IK . EQ . 3 ) Z ( I ) *250 
I F  f  IK.  E0.3)TU)«.  83333 
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IF  (  IK.EQ.4)X(n>r  CR0SSM-U*.5'  +  (CROSS  '1  '*.833  33  ' 

IF  ( IK.EQ.4)Y(  I  )«(RANGE(I-1)*.5  WRANGEf  I)  *.83333  ^ 

IF( IK.E0.4 )Z( I )«400. 

IF(IK.EQ.4)T(I )*.5+. 83333 
I F (IK . LE . 4  )GO  TO  21 
33  CONTINUE 

IFflK.GT. 12)60  TO  44 

X ( I )«(CROSS(I-l)*.83333)*(CROS$( I  )*  .83333) 

Y (I )- (RANGE (  1-1)*. 833 33) ♦( RANGE ( I)*.  83333  ) 

Z ( 1 ) «  500 . 

T(I)«1. 66667 
GO  TO  21 
44  CONTINUE 

IF ( IK . EQ . 13 ) X ( I )•( CROSS ( I- 1)*. 83333)* (CROSS (I )*1 .66667) 
IF ( IK . EQ . 13 ) Y ( I ) -(RANGE! I - 1  )*. 83333 )♦( RANGE ( P *1 . 66667  ) 
IF(IK.EQ.13)Z(I)»750. 

I F { IK . EQ . 1 3 ) T ( I )■. 83333+1. 66667 
IFdK.EO. 13)60  TO  21 

X( I )« (CROSS ( I -1)*1. 66667 )f(CROSS (I)*1.66667) 

Y ( I )»(RANGE( I-  1  1  *1.66667  ) ♦( RANGE (  I  ) *1.66667  ) 

Z ( I ) » 1000 . 

T( I ) «3 . 33334 
21  CONTINUE 

C 
C 

C  THE  PROGRAM  USES  THE  TIME  AND  POSITION  OF  THE  FIRST 
C  STATION  READ  IN  AS  THE  ORIGIN,  THE  FOLLOWING  COMPUTES  THE 
C  INITIAL  OFFSET  FOR  EACH  SUCCESSIVE  STATION 
C 
C 

IF( ( IK.EQ.l) .AND. (K . GT  ,  1 } )  X (  I  )»CXfK)-CX(l) 
IF((IK.EQ.1).AND.(K.GT,1))Y(I)»CY(K)-CY(1) 

IF( (IK.EQ.l ).AND.(K.GT.l))Z(l)«CZ(X)-CZ( 11 
IF  (( IK.EQ.l).  AND.  (K.6T.1))DELT«Z(  I  W1./300.) 

I F( (IK.EQ.l) .AND. (K.GT.l) )TH» ( THR (K ) -THR ( 1)1*60. 

IF(  (IK.EQ.l) .ANO.(K.GT. 1) ) TM» ( TMI N (K ) -TMI N ( 1 )  ) 

IF ( ( IK.EQ.l) .AND. (K.GT.l) )T( I 1-TH+TM-DELT 
IF(IK.GT.1)X(I )-X(I)*X(I-l) 

IF(IK.GT.1)Y(I)»Y(I)+Y(I-1) 

IFflK.GT. 1)Z(I)»Z(I)+Z(I-1) 

IF(IK.6T.1)T(!)«T(I)*T(I-1) 

IZ»IZ+1 

20  CONTINUE 
RETURN 
END 

0MAP.IL  , PREDICT/A 
LIB  NR-A*RBL I B  . 


sample  data  deck  using  four  stations  two  times 
and  ten  layers. 


4  2  10 


503109. 

189735. 

4033  . 

10. 

00. 

557402. 

189530. 

4198. 

10. 

00. 

472572. 

215268. 

3999. 

10. 

00. 

543736. 

140375. 

4097. 

10. 

00. 

503109. 

189735. 

4033. 

12. 

00. 

557402. 

189530. 

4198. 

12. 

00. 

472572. 

215268. 

3999. 

12. 

00. 

543736, 

140375. 

4097. 

12. 

00. 

11151000  TSX 
10 

11151000  0R0 
10 

11151000  SMR 
10 

11151000  MCG 
10 

11151200  TSX 
10 

11151200  ORO 
10 

11151200  SMR 
10 

11151200  MCG 
10 

010005 
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SAMPLE  OUTPUT 


@XQT  PREDICT/A 
@ADD  FL. DAT A/ NEW 


503109. 

189735. 

4033. 

10. 

0. 

557402. 

189530. 

4198. 

10. 

0. 

472572. 

215268. 

3999. 

10. 

0. 

543736. 

140375. 

4097. 

10. 

0. 

503109. 

189735. 

4033. 

12. 

0. 

557402. 

189530. 

4198. 

12. 

0. 

472572. 

215268. 

3999. 

12. 

0. 

543736. 

140375. 

4097. 

12. 

0. 

1 

4730.000 

6.000 

288.3 

875 

2 

3660.000 

2.000 

290.7 

865 

3 

4200.000 

4.000 

288.6 

840 

4 

4650.000 

13.000 

284.8 

810 

5 

5020.005 

16.000 

285.0 

745 

6 

4724.610 

38.209 

283.2 

708 

7 

5329.776 

43.313 

279.8 

666 

8 

4670.980 

41.689 

276.0 

626 

9 

4802.574 

60.206 

272.2 

589 

10 

4816.960 

57.369 

268.2 

554 

46 


